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ABSTRACT 

In  this  report,  we  have  designed  an  essentially  non-oscillatory  reconstruction  for  functions 
defined  on  finite-element  type  meshes.  Two  related  problems  are  studied  :  the  interpolation 
of  possibly  unsmooth  multivariate  functions  on  arbitrary  meshes  and  the  reconstruction 
of  a  function  from  its  average  in  the  control  volumes  surrounding  the  nodes  of  the  mesh. 
Concerning  the  first  problem,  we  have  studied  the  behaviour  of  the  highest  coefficients  of  the 
Lagrange  interpolation  function  which  may  admit  discontinuities  of  locally  regular  curves. 
This  enables  us  to  choose  the  best  stencil  for  the  interpolation.  The  choice  of  the  smallest 
possible  number  of  stencils  is  addressed.  Concerning  the  reconstruction  problem,  because  of 
the  very  nature  of  the  mesh,  the  only  method  that  may  work  is  the  so  called  reconstruction 
via  deconvolution  method.  Unfortunately,  it  is  well  suited  only  for  regular  meshes  as  we 
show,  but  we  also  show  how  to  overcome  this  difficulty.  The  global  method  has  the  expected 
order  of  accuracy  but  is  conservative  up  to  a  high  order  quadrature  formula  only. 

Some  numerical  examples  are  given  which  demonstrate  the  efficiency  of  the  method.  _ 
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1  Introduction 


During  the  past  few  years,  a  growing  interest  has  emerged  for  building  high  order  accurate 
schemes  (i.e  of  order  greater  that  2)  for  compressible  flows  simulations.  It  is  well  known  that 
even  for  smooth  initial  conditions,  these  flows  may  develop  discontinuities  that  make  linear 
schemes  useless. 

At  the  beginning  of  the  80’s,  the  class  of  Total  Variation  Diminishing  schemes  appeared 
and  they  have  been  successfully  and  widely  used  with  many  types  of  meshes  (see  for  example, 

[1]  for  a  review  and,  among  many  others,  [2]  for  simulations  on  finite  element  type  meshes). 
Nevertheless,  one  of  their  main  weaknesses  is  that  the  order  of  accuracy  falls  to  first  order 
in  regions  of  discontinuity  and  at  extrema,  leading  to  excessive  numerical  dissipation. 

Various  methods  have  been  proposed  to  overcome  this  difficulty  (adaptation  of  the  mesh 
for  example)  but  one  promising  way  may  also  be  the  class  of  the  Essentially  Non-Oscillatory 
schemes  (E.N.O.  for  short)  introduced  by  Marten,  Osher  and  others  [3,  4,  5,  6,  7]. 

The  basic  idea  of  E.N.O  schemes  is  to  use  a  Lagrange  type  interpolation  with  an  adapted 
stencil  :  when  a  discontinuity  is  detected,  the  procedure  looks  for  the  region  around  this 
discontinuity  where  the  function  is  the  smoothest.  Then  a  reconstruction  technique  may  be 
applied  which  enables  approximation  of  the  function  to  any  desired  order  of  accuracy  from 
its  average  in  control  volumes  surrounding  the  mesh  points.  The  approximation  is  done  so 
that  it  is  conservative. 

Some  attempts  have  been  made  to  extend  these  ideas  to  multidimensional  flows  (see  for 
example  [8]),  but  only  for  structured  meshes. 

In  this  report,  we  intend  to  study  the  problem  of  the  reconstruction,  up  to  any  order 
of  accuracy,  of  a  given  function  given  either  by  its  value  at  the  nodes  of  a  triangulation  or 
by  its  averages  on  control  volumes  defined  around  these  nodes  so  that,  in  the  second  caise, 
the  reconstruction  is  conservative.  This  latter  problem  has  already  been  studied,  for  smooth 
functions  only,  by  Barth  et  al.  [9]  but  their  method  does  not  appear  to  generalize  easily  to 
unsmooth  functions. 

The  outline  of  this  report  is  as  follows.  In  the  first  part,  we  give  some  basic  facts 
about  Lagrange  interpolation  in  several  dimensions.  In  particular,  we  study  the  problem  of 
the  localization  of  regions  of  smoothness  from  the  Lagrange  interpolation  coeflRcients  that 
generalize  those  known  in  one  dimension.  These  results  seem  to  be  new.  Then,  we  propose 
an  algorithm  for  E.N.O.  interpolation  and  we  try  to  give  some  indications  for  selection  of  the 
smallest  family  of  the  possible  stencils  for  second  and  third  order  approximation,  the  only 
ones  considered  in  this  report.  We  also  propose  an  adaptation  of  the  so-called  reconstruction 
via  deconvolution  procedure  that  was  originally  built  for  regular  meshes,  and  indicate  why, 

or 
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in  general,  the  conservation  property  must  be  lost  to  ensure  a  high  order  of  approximation. 
Some  numerical  tests  indicate  the  performance  of  this  method. 


1.1  Notations 

•  IR„[X,  F]  :  finite  dimensional  vector  space  of  iwo  variable  polynomials  over  IR, 

•  N{n)  =  ^  ■  dimension  of  IR„[X,  F], 

•  admissible  stencil  for  solving  the  Lagrange  nroblem  in  1R„[X,  F],  see  section  2.1.2, 

•  ||f7||  is  the  euclidian  norm  of  U, 

•  £)”u  :  n-th  derivative  of  u. 


2  Lagrange  interpolation  on  arbitrary  sets  of  IR^ 

2.1  Sets  of  admissible  points 

2.1.1  The  polynomials  in  IR^ 


We  will  denote  by  IRfX,  F]  the  vector  space  of  the  polynomials  of  two  variables  (X  and  F) 
with  coeflRcients  in  IR.  An  element  of  IRfA",  F]  may  be  described  by  its  (finite)  expansion  in 
terms  of  powers  of  X  and  F  ; 


/=1  i+j=l  i,j>0 


(1) 


The  highest  integer  such  that  at  least  one  of  the  coefficients  of  the  monomials  is  non 

zero  is  called  the  tota!  degree  of  P. 

If  (xo,j/o)  is  a  point  of  IR^,  another  expansion  of  P  may  be  written  in  terms  of  the 
monomials  (X  —  xo)’(F  —  yoY  with  the  help  of  the  Taylor  formula.  The  total  degree  of  P 
does  not  depend  on  the  point  (xo,j/o)- 

In  the  sequel,  we  will  denote  by  IRn[A’,F]  the  (finite)  vector  space  of  the  polynomials 
of  IR[A’,  F]  with  total  degree  less  or  equal  to  n.  This  vector  space  has  dimension  N{n)  = 
(n  +  l)^(n  +  2)^  ^  basis  of  which  is  the  set  of  monomials  {X  —  Xo)‘(F  —  yoY  of  total  degree 
i  +  j  less  or  equal  to  n. 

Let  us  now  describe  another  interesting  basis  of  IRnfA",  F].  Consider  {A,  5,  C)  a  triangle  of 
IR^  and  let  us  denote  by  Aa,  Ag,  Ac  the  barycentric  coordinates  of  the  three  points  (A,  B,  C) 


2 


defined,  for  any  point  M,  by  : 


M  —  AjX  a  +  Ab  B  +  Ac  C 
Aa  +  Ab  +  Ac  =  1 


(2) 


It  is  easy  to  se  that,  for  any  pair  of  points,  say  A  and  B,  the  set  A;^ 
basis  of  IRnlA",  Y]. 


K+j<n 


is  also  a 


2.1.2  The  Lagrange  interpolation  problem  in  IR^ 

The  Lagrange  interpolation  problem  may  be  formulated  as  follows  : 

Given  N  and  n,  two  integers,  a  family  of  N  points  in  1R^  and 

N  real  values  (wi)i<j<;yf,  find  an  element  P  o/IR„[X,  y]  such  that  for  any  point 
Ai  of  one  has  P(Ai)  =  Ui. 

In  the  sequel,  we  will  often  make  no  distinction  between  an  element  of  Ai,  and  its 
coordinates  in  a  suitable  frame,  {xi,yi). 

For  this  problem  to  have  a  solution,  two  conditions  must  be  fulfilled  : 

1.  one  must  have  N  = 

2.  the  following  generalized  Van  der  Monde  determinant  must  be  non  zero  ; 


A^cn)  =  det  [x)  yf] 

i+j<n 

1  xi  yi 

•••  a:?  X?  Vi  ' 

•••X? 

ixi,yi)  e 

1  Xff  yN 

•  •  •  x]y  ^y^v  • 

••x]V 

We  will  say  that  the  set  is  admissible  if  A5(n)  7^  0.  If  a  set  is  admissible,  there 
exist  ^  coefficients  ,  (a,  j),  such  that  the  solution  of  the  Lagrange  problem  is  : 

^  =  E  E  (4) 

l=l  i+j=l  i,j>0 

The  problem  of  characterizing  the  admissible  sets  has  been  widely  studied,  see  [10]  for 
example  and  the  references  therein. 

Remarks  : 

1.  The  condition  (3)  has  been  given  for  the  basis  X'Y^  of  IRnfA",  Vj.  A  similar  and 
equivalent  condition  could  have  been  given  for  the  two  other  bases  we  have  mentioned . 
The  formula  (4),  provided  that  the  monomials  X'Y->  are  replaced  by  the  elements  of 
the  new  basis,  is  also  true. 
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2. 


If  card  <5^*^  =  3,  this  condition  is  nothing  more  than  the  one  which  says  that  the  three 
points  must  not  be  aligned. 


3. 

4. 


If  we  were  in  IR,,  this  determinant  would  be  the  classical  van  der  Monde  determinant. 

The  set  of  ^  ^'^-uplets  where  the  condition  (3)  is  not  fulfilled  is  an  algebraic 

curve  of  and  consequently  a  closed  subset  of  measure  zero  in  lR7V(n) 


In  the  next  section,  we  address  the  question  of  the  practical  calculation  of  the  coefficients 


2.1.3  Determination  of  the  Lagrange  expansion 

In  this  section,  we  use  the  monomials  X'Y^  for  expanding  polynomials,  but  any  other  basis 
would  be  suitable  and  the  results  are  immediately  transferable. 

The  coefficients  of  the  Lagrange  expansion  are  the  solution  of  the  linear  N{vi)  x  N{n) 
system  : 

ui  =  'f2  "«■  yi  ^  (5) 

l=l  i+j=l  i,j>0 

Since  condition  (3)  is  true,  the  Cramer  formula  applied  to  (5)  gives  the  answer. 

Several  authors  have  tried  to  generalize  the  Newton  formula  that  make  the  Lagrange 
interpolation  efficient  from  a  numerical  point  of  view,  and  a  very  general  answer  has  been 
given  by  Muhlbach  [11,  12]. 

In  these  papers,  he  addresses  the  problem  of  the  Lagrange  interpolation  by  a  set  of 
functions  on  a  set  of  points  He  calls  the  set  (/i)  a  Cebysev- system  if  given 

any  function  /,  for  any  pair  of  subsets  of  /,  L  and  M,  having  the  same  (finite)  number  of 
elements,  there  exist  real  numbers  Qj  such  that  : 


Ui  =  ^  ajfj{Ai),  for  all  Ai  G  L. 

jeM 


(6) 


For  the  sake  of  clarity,  we  may  assume  that  L  =  M  =  {!,•  •  •  N}.  He  uses  the  notation  : 

/ 


A\  -  ■  ■  Ak 


for  denoting  the  coefficients  of  fk  in  the  development  (6).  Then,  in  [12],  he  shows  that  if  one 
has  a  Cebysev  system  (theorem  4.1  pp.  106)  : 


Ai  ■  ■  ■  An 


f 


/l  •  •  •  /n-1 

A2--'  An 

/]- 

/l  • 
Ai  • 

•  /n-1 

•  •  -^n-l 

/I 

/l  •  •  •  /n-1 

A2---  An 

/I- 

■/i-- 

•  /n-1 
■  •'^n— 1 

/■ 

(7) 
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This  expression  is  a  direct  generalization  of  the  classical  Newton  formula.  Let  us  make 
several  comments  on  this  formula  when  applied  to  our  problem  : 

1.  If  one  adopts  the  lexicographic  ordering  1,X,Y ,X^,XY,Y^,- •  ■ 

•  •  •,XY’^~^ ,Y"',  the  previous  formula  (7)  must  be  applied  n  +  1  times  to  go  from  a 
total  degree  n  to  a  total  degre  n  +  1.  One  must  also  store  quite  a  lot  of  terms  like 

/i  •  ■  •  /fc  f 

Ai‘--Ak  ^ 

to  build  the  divided  difference  table.  For  example,  to  go  from  degree  one  to  degree 
two,  one  must  evaluate  and  store  C|  approximations  of  gradients  by  means  of  approx¬ 
imations  on  triangles,  combine  them  to  obtain  approximations  on  sets  of  four  points 
(C|  sets),  of  five  points  (Cg  sets)  and  of  six  points  (one  set).  Moreover,  to  go  from 
approximation  on  k  points  to  fc  -|-  1  points,  (k  +  1)  x  (k  +  1)  determinants  must  be 
evaluated. 

2.  From  a  numerical  point  of  view,  the  basis  X'Y^  or  {X  —  a)‘(i^  —  by  are  not  well  suited 
to  calculations.  This  can  easily  be  seen  since  for  any  pair  {a,b), 

H*'']  i+j<n  =[(»'-“)’ (W-Vl I  j  +  j<„  • 

If  (a,  6)  is  any  point  of  and  if  /i  =  j,j)g5(n)(lx(  —  a|,  |r//  — 6]),  then  Hadamard’s 

inequality  shows  that  : 

|A5(„)1  < 

where  K(n)  =  1  -f  Y,i=i,n  ^  =  0{rA)  so  that  one  reaches  very  quickly 

machine  zero  though  the  linear  system  may  be  well  conditioned. 

An  alternative  to  this  last  point  is  to  use  local  coordinates  such  as  the  barycentric  ones.  In  the 
E.N.O.  method  we  will  develop  in  a  further  section,  for  each  point,  the  natural  barycentric 
coordinates  are  not  known  a  priori  so  that  the  work  has  to  be  repeated  at  each  interpolation 
call.  If  this  is  included  in  an  iterative  algorithm,  the  cost  (and  the  storage)  seems  to  be  much 
too  important  at  least  for  the  cases  we  have  considered  in  this  report.  For  all  these  reasons, 
we  have  preferred  to  use  classical  inversion  techniques  for  linear  systems. 

2.2  A  recurrence  formula 

In  this  section,  we  wish  to  show  another  recurrence  formula  that  enables  us  to  obtain  the 
coefficients  of  the  expansion  of  total  degree  n  from  those  of  the  expansion  of  total  degre  less 


than  n  in  only  one  step.  This  recurrence  formula  may  be  viewed  as  another  version  of  that 
given  in  [12],  theorem  3.1,  page  400  and  will  be  useful  in  section  2.3. 

Let  us  begin  with  some  notations.  Let  P\,  •  •  •,  Ps(n)  be  a  basis  of  1R„[X,  T]  such  that 
Pi,  ■  •  •)  Pn(p),  P  <  n,  is  a  basis  of  IElp[X, K].  The  three  bases  we  have  considered  in  section 
2.1.2  are  of  that  kind.  Let  Ai,  •  •  •,  A^^n)  be  admissible  points.  We  set  : 

Rl  =  (Pl(Ar,)  •  •  •  PN(n)(AL)) 

so  that  the  solution  of  the  Lagrange  problem  (where  the  u,’s  are  given), 

u,  =  ^  ctjPj{Ai)  for  all  Ai  £  L, 

l<j<N(n) 

may  be  seen  as  the  solution  of  the  linear  system 

M  (ai---  a^^n))  =U  -  (ui  •  •  •  UN(n^  ,  (8) 

where  the  Lth  row  of  M.  is  Ri.  The  solutions  of  system  (8)  are  : 


„  _  det{Ri  •••  Rl'-'  RN(n)) 
^  det{Ri  •  Rl  '  RN(n))  ’ 


where  Rl  =  U  denotes  the  L-th  column. 


(9) 


Lemma  2.1  Let  (A,)i<,<yv(n)  be  an  admissible  set  in  which  any  of  its  N{p)^p  <  n  subsets 
is  admissible.  Then,  for  a  given  p  <  n,  let  1  =  {ii,  •  •  •  ?;v(p)}  a.nd  J  he  ordered  sets  such  that 
/U  J=  {!,•  -  ■  N{n)}.  If  A  =  (a,-  j)  is  a  N{n)  x  N{n)  matrix,  we  set 

det{A)i  =  det{ai  j)  ^  <  ^-  <  ;v(p)  ’ 

J  e  I 

and  det{A)j  =  det(G. ,)  +  i<i<  N{n)  ’  ’ 

Let  us  also  denote  by  the  coefficient  of  Pv  of  the  Lagrange  problem  of  degree  p  for  nodes 
in  I  (for  degree  n,  we  omit  the  subscrjipt  I). 

Then  for  any  U  <  N{p),  we  have 

=  lll,caTd{I)=N(p)  ^  I 

where  a{I)  =  1  +  p{p  +  l)/2  +  (^O)i  Rl.'  appears  a  first  time  at  the  L'-th  row  of 

det{Ri  ■  •  •  Ry  •  •  ■  Rn')i  and  a  second  time  at  the  L—N{p)-th  row  of  det{RN ■  •  •  Ri,,  •  •  •  Rn)j. 


-  r-i 


det{Ri 


Rni)i  det{RN{p)+i 
det{Ri  •  •  •  Ri^) 
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Proof  :  By  switching  columns  L  and  L*  in  (9),  one  gets  : 


det{Ri  •  •  •  Rl  •  •  •  Ri!  •  •  •  RN(n)) 
det(Ri  •••/?£,•••  RN{n)) 


Then,  a  direct  application  of  the  generalized  Lagrange  formula  (see  [13],  pp.  19-22)  to  the 
previous  expression  gives 


det{Ri,  ■  Rl  -  •  •  Rv  -  -  •  RN(n))  = 

12l,card(I)=N(p){~^Y^^^  dtt{R\  •  ■  •  Rl  •  •  •  Rn{p))i  del(i?jv(p)+l  •  •  •  Rl'  •  -  ■  RN(n))j 


Since  any  N{p)  subset  is  admissible,  one  has 

(p)  _  det{Ri  Rl---  Rn(p))i 
^  ^  dei{Ri  ■  •  •  Rl’  -  •  •  -Rjv(n))/  ’ 

and  the  results  follows  immediately.  • 


Lemma  2.2  With  the  assumptions  of  lemma  2.1,  then,  forp  <  n,  if  N{p)  <  L  <  N{n)  and 
V  <  N(p), 

E  At"'=o. 

I,card(l)=N(p) 

Proof  ;  Apply  the  Lagrange  formula  to 

det{Ri  •  •  -  Rli  -  •  -  Rn(p)Rn{p)+i  •  •  •  RL’RN(n))  =  0, 

and  interpret  the  coetBcients  in  terms  of  q’s,  all  equal  to  1,  and  A’s.  The  lemma  2.1  gives 
the  result  • 


2.3  Approximation  of  smooth  and  unsmooth  function  by  poly¬ 
nomials 

The  problem  of  interest  in  this  section  is  the  following  :  Let  u  be  a  real  function  defined 
on  an  open  subset  fl  of  We  assume  that  u  is  n  times  continuously  differentiable  on  fl 
except  perhaps  on  a  subset  of  consisting  of  a  finite  collection  of  locally  curves.  Let  now 
T  be  a  mesh.  For  each  point  of  T,  we  consider  a  Lagrange  interpolation  of  u.  Is  it  possible 
to  localize  the  regions  of  smoothness  of  u  from  the  coefficients  of  the  Lagrange  interpolation 
of  u  ?  The  answer  is  yes,  at  least  for  second  and  third  order  approximations  if  additional 
iissumptions  are  made  on  the  mesh.  These  assumptions  guarantee  that  one  can  solve  the 
Lagrange  problem  for  any  order  from  1  to  n,  and  may  be  seen  as  a  very  natural  generalization 
of  classical  conditions  used  in  the  finite  elements  theory  [14]. 

For  functions  defined  on  IR,  one  knows  that  the  divided  differences  of  u  satisfy  : 


7 


•  If  u  is  smooth  on  an  interval  I  containing  xi,  •  •  •  x„,  then  there  exists  ^  €  /  such  that 

r  11- 

n! 

•  if  has  a  jump  on  /,  one  ha^ 

[zi,  X2,  •  •  •  ,  Z„|u]  =  0([u('')](Zn  -  Xi)*""). 

In  this  section,  we  intend  to  generalize  these  relations,  and  in  particular,  the  second  one 
since  this  problem  seems  (surprisingly)  not  to  have  been  studied  yet.  The  proof  appears  to 
be  technical,  and  we  have  not  been  able  to  prove  it  for  any  total  degree.  The  proof  is  divided 
into  two  parts.  In  the  first  part,  we  study  the  case  of  a  stencil  of  N{n)  points  where  u 
admits  two  values,  0  and  1.  We  show  here  that  the  polynomial  of  degree  n  that  interpolates 
u  is  exactly  of  total  degree  n.  Then,  we  define  a  condition  on  the  stencils  that  appears  to  be 
a  generalization  of  the  one  that  says  that  triangles  must  not  have  too  small  angles  to  ensure 
a  uniform  error  bound  for  classical  finite  elements  [14].  Then,  using  Lemmas  2.1  and  2.2, 
we  obtain  our  result.  Let  us  begin  with  the  case  of  a  stencil  in  which  the  convex  hull  u  is 
smooth. 

2.3.1  Case  of  a  “smooth”  stencil 

This  problem  has  been  studied  by  for  example  Ciarlet  and  Raviart  in  [15].  Let  us  recall  one 
of  their  main  results  : 

Theorem  2.3  Let  S  =  be  an  admissible  (for  degree  k)  set  of  points  o/IR”,  and  let 

h  and  p  be  respectively  the  diameter  of  E  and  the  suprenum  of  the  diameters  of  the  spheres 
contained  in  the  convex  envelop  K  ofE.  Let  u  be  a  function  that  admits  everywhere  in  K  a 
k  +  1th  derivative  with 

Mk+i  =  sup{||i9*‘''^u(z)||;x  €  K}  <  +00. 

If  P  is  the  unique  interpolating  polynomial  of  degree  <  k  of  u,  we  have  for  any  integer  in 
with  0  <  m  <  k, 

sup{||£>"‘«(x)  -  D^Pix)\\-,xe  K)  < 

P 

for  some  constants 

C  =  C(n,  k,  m,  S). 

Moreover,  if  E'  is  obtained  from  E  by  an  affine  transformation,  then  C{n,  k,m,E)  = 

C(n,  k,m,  S'). 

From  this  inequality,  one  sees  that  the  “flatter”  E  is,  the  poorer  the  estimation  is.  A 
direct  application  of  this  theorem  gives  a  generalization  of  our  first  statement. 


8 


2.3.2  Case  of  an  “unsmooth”  stencil 


Study  of  a  simplified  problem 


Let  us  consider  an  admissible  stencil  of  cardinality  Sq,S\  two  non¬ 

empty  subsets  of  5^"^  having  empty  intersection,  the  union  of  which  is  Let  us  consider 
a  polynomial  P  of  total  degree  n  such  that  for  all  points  of  5o,  P  has  value  0  and  for  those 
of  «5i,  P  has  value  1.  Then  we  conjecture  that  : 

Conjecture  2.1  If  is  admissible,  then  the  total  degree  of  P  is  exactly  n. 

We  have  not  found  in  the  literature  any  general  proof  of  the  statement,  but  we  have  the 
following  lemma  : 

T  _ n  A  rr _ .r _ j  1-1  +  +  ^ 


Lemma  2.4  If  any  subset  of  cardinality 
the  conjecture  2.1  is  true  for  n  =  1,2 


k  <  n,  of  is  admissible,  then 


Proof  ;  The  proof  will  be  given  for  degree  1,2.  We  also  indicate  the  difficulties  for  higher 
degree. 

•  Degree  one.  The  stencil  is  made  of  three  points  A,  B,  C  that  form  a  triangle.  P  is 
either  of  type  or  1  —  and  is  of  degree  exactly  one. 

•  Degree  >  1.  Let  us  assume  that  P  is  at  most  of  degree  n  —  1.  Set  N  =  card{So)  and 
A/  =  <5i.  We  have  N  +  M  =  N{n).  One  may  assume  that  N  <  M  by  changing  P  into 
1  —  P.  So,  2N  <  N{n).  In  the  following  table,  we  give  the  maximum  values  of  N  up 
to  degree  6  : 


degree 

2 

6 

3 

3 

10 

5 

4 

15 

7-8 

5 

21 

10-11 

6 

28 

14 

From  this  table,  one  can  see  that  there  is  always,  for  degree  2,  at  least  3  points  that 
have  the  same  value.  Since  these  three  points  are  admissible  for  degree  one  and  since 
by  assumption,  P  is  either  a  constant  or  of  degree  1,  we  see  that  it  must  be  a  constant 
which  is  absurd  since  it  takes  two  different  values.  The  same  argument  applied  to 
degrees  n=3,4,5  shows  that  P  must  be  of  degree  exactly  n  —  1  because  we  always  have 

(n  -f-  1)  (n  -f  2)  ^  n  (u  -I-  1) 

4  “  2  ’ 
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but  fails  for  degrees  greater  or  equal  to  6  (because  the  previous  inequality  does  not 
hold  if  n  >  5).« 

With  this  in  hand,  we  get  the  following  result,  if  Ri  j  is  the  following  vector  : 

j  =  ((3=1  -  a:o)‘(3/i  -  yoY  •  ■  •  (xjv  -  XoYivN  -  yoY)^ . 

Lemma  2.5  Let  (xoij/o)  be  any  point  of  the  convex  hull  E  of  Assume  that  conjecture 
2.1  is  true.  Let  P  be  a  polynomial  that  is  1  on  Si  and  0  on  So  where  both  sets  satisfy  the 
assumptions  of  lemma  U.j.  If  h  =  max{{\xi  -  xo\,\yi  -  yQ\),[xi,yi)  e  5^"^},  and  t/ 5^"! 
satisfies,  for  some  a  >  0, 

then  there  exist  two  constants  Ci{n,a)  and  C2{n)  such  that  the  coefficients  of  the  Taylor 
expansion  of  P  around  (xo,yo)  : 

P  =  £  a,  j(X-  x„)‘(y  -  noY, 

«+j<n 

satisfy 

C2(n)k-->  Y.  \^il\>Ciin,a)h-\  (12) 

«+j=n 

Remark  : 


1.  The  set  of  points  we  define  is  clearly  not  empty  because,  on  the  convex  hull  of  5^”^ 
(which  is  compact),  the  function  defined  by  the  right  hand  side  of  (11)  is  continuous. 

2.  If5W  were  a  triangle,  for  example,  the  minimum  of  that  function  would  be  the  mini¬ 
mum  of  the  absolute  value  of  the  sines  of  its  angles. 

Proof  :  We  adopt  the  lexicographic  ordering  of  monomials.  The  set  of  admissible  stencils 
satisfying  condition  (11)  and  j/ij  <  C  is  a  compact  subset  C  of  Let  us  consider  the 

real  functions  defined  on  C  by  : 

4>i  •  •  • )  An)  =  Yf  ■  ■  -  Ui  j  -  ■  •  i?On)|  7 


I  /  A  /I  \  detl^Roo  ■  ■  ’  Vij  •  •  •  Ron) 

'  det(R,-R,)-  • 

in  which  the  vector  Uij  stands  at  the  “ij”-th  column  and  has  the  value 


L-,  =  (0,-.-,l,-..0)' 
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Let  also  V<  j  be  P{Aij)Uij  where  the  “1”  is  at  the  L-th  position  (we  refer  to  the 

lexicographic  order).  Let  us  note  that  P{Aij)  is  zero  or  one,  and  its  value  depends  only  on 
ij  and  not  5^"^ 

It  is  clear  that  ipij  is  the  coefficient  of  {X  —  xoy{Y  —  yoY  in  the  Taylor  expansion  of  P. 
It  is  also  clear  that 

\det[RQ  0  •  •  •  i  •  •  •  n]l 


by  the  triangle  inequality.  So, 


=  Y.  I“mI>  E 

i+j=n  t+j=n 


dc  t  (  Ftoo  *  *  *  9  *  *  *  Ron ) 

4>i  j 


The  left  hand  side  of  this  inequality  is  a  continuous  function  on  C  and  hence  reaches  its 
minimum.  This  minimum  cannot  be  zero  because  that  would  mean  that  all  of  the  a,  j’s  are 
zero  which  is  absurd  by  assumption,  and  from  lemma  2.2. 

Now  let  us  turn  to  the  second  inequality.  We  have 

A"«(A,..  -.A«)  =  A”  ■£,  E  i*il- 

•+j=n  »+j=n 

The  left  hand  side  is  also  a  continuous  function  on  C  and  is  bounded  above.  Clearly,  the 
latter  constant  does  not  depend  on  o;.  • 

Corollar  2.6  With  the  assumptions  of  lemmas  2.4  and  2.5,  let  n  and  p  be  integers  satisfying 
p  <  n.  Choose  L  and  L'  such  that  N{p  —  1)  <  L'  <  N{p)  and  N{n  —  1)  <  L  <  N{n).  Then 
there  exist  two  constants  C\  and  C2  such  that 

C2(n)/i-"+'’  >\^^'>  (7i(n,a)/i-"+P, 

for  any  subset  of  cardinality  p. 

Proof  :  We  apply  the  definition  of  Af  and  the  same  techniques  used  in  the  previous 
lemma  • 

Now,  we  may  state  our  main  result  : 

Theorem  2.7  Let  be  a  stencil  satisfying  the  assumption  of  Lemmas  2.4  and  2.5  for 
n  =  1,2.  Let  u  a  real  function  defined  on  an  open  subset  ofU  in  IR^  being  except  perhaps 
on  a  locally  curve  where  its  n-th  order  derivative  may  have  a  jump  [D"u]  such  that  the 
intersection  I  of  that  curve  and  the  convex  hull  of  is  not  empty.  Then  there  exits  a 
constant  C(n,a)  >  0  such  that  the  coeffirients  in  the  Taylor  expansion  towards  a  point  ofJ 
satisfy  : 

X;  \ai  j\>C{n,a)^-^f^.  (13) 

i+J=n  ^ 
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Proof  :  Let  us  assume  that  u  admits  p-th  continuous  derivatives  but  its  p+  1-th  ones  have 
a  jump  on  a  locally  curve  C.  Let  be  the  coefficients  of  a  Lagrange  interpolation  of 
degree  p  on  a  suitable  subset  of  5^”^.  The  recurrence  formula  of  Lemma  2.1  gives  : 

E 

/,card/=N' 

Let  S+  and  <S_  be  the  subset  of  defined  by 


a  ^  S+{  resp.  5_)  iff  Diju{a)  — >  £)+(  resp.  D-^) 


These  two  sets  exist  because  of  the  topological  nature  of  the  curve  C  (see  Figure  1). 

Let  e  >  0  and  xq  be  a  point  of  J.  There  exists  Tj  >  0  such  that  for  max{\ax  —  xqI,  \ax  — 
^ol)  <  we  have  the  following  development  of  u  for  points  of  5+  : 

“(«)  =  E  H  Diju{xo){ax  -  xoYiUy  -  yoY  +  ^  +  o(l))  (a^  -  XoYiay  -  J/o)^ 

/=0  i+j=l  t-l-j=p+l 


where  Ox  and  Oy  are  the  x~  and  y-  components  of  a  and  Dfj  is  the  limit  of  the  i  j  partial 
derivative  of  u  when  x  xq  on  the  right  and  |o(l)|  <  e  .  The  same  type  of  development  is 
also  true  for  points  of  5_  ;  Dfj  is  replaced  by  D"  .  The  properties  of  the  determinants  allow 
us  to  get  for  any  subset  of  cardinality  N{p)  that,  by  changing  the  ordering  of  the  points  if 
necessary, 


xoYial-yoy 

Dt,{a7-Xoy{a^-yoy 

D-^ia^+^-xoYia^^^-yoy  ■  ■  ■  Bop 

D-ya^^»^-xoYia':!^’’^-yoy _ 

det{Roo  ■  ■  ■  i^p) 
o(l)(ai  -  a;o)’(aJ  -  yoY 

-  a:o)'«  -  yoY 

/2oo  •••  o(l)(ar‘-^o)’(a;r^'-yoF  ' '  ’ 

o(l)(af  P)  -  -  yoY 

det{Roo  •  •  •  Rop) 


(14) 


where  the  index  L'  stands  for  {i,j)  in  the  lexicographic  order.  For  the  sake  of  sim.plicity,  let 
us  denote  by  p/  and  i//  the  coefficients  obtained  from  the  first  part  of  the  right  hand  side  of 
equation  14  by  replacing  £)+  by  1  and  D~j  by  0  for  p/  and  vice  versa  for  vi  so  that 


“M  =  [Dtj  +»(!))/•/+  [DJi  +  0(1)) 
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We  have  to  notice  that  the  sum  of  ni  and  vj  is  one,  and  that  if  all  the  points  of  /  are  on  the 
same  side  of  C,  then  either  /x/  or  ui  is  zero.  Then,  using  lemma  2.2  and  the  above  remarks, 
we  have,  if  [Dij]  =  Dfj  —  D~j  : 

=  [Aj]  X)  IT  ^ 

card{I)=p  card(I)=p  card(I)=p 


SO 


that 


\[Dij]\ 


L  L' 


caTd(I)=p 


<M  +  £  E  (I/*;!  +  k/l)  |Af ‘'I. 


card(I)=p 

We  can  also  get  a  similar  inequality  by  replacing  m  by  i/j  if  the  points  of  I  are  not  on  the 
same  side  of  C  so  that  : 

A'-”  <  A--"  {KAill  |E„..(;).pWAf 


+  l[Aill  |E„,j(;)=,<';Af '■'D 

+  2  ftf-latl  +  (b/l  +  k/l)  |Af  ‘-'I 

(15) 

When  all  the  points  of  I  are  on  the  same  side  of  C  the  factor  2  is  replaced  by  1 .  Because  of 
inequality  (11)  and  from  Hadamard’s  inequality,  one  has  : 


l/^/l  +  Wi\  < 


-  yoY^  +  ./5I  (oi  -  -  yoy^ 

Y  a  S+  Y  ® 


a\\Ru 


This  latter  expression  is  bounded  above  by  a  constant  C  because  L'  corresponds  to  (*,j) 
in  the  lexicographic  order,  without  any  additional  conditions  on  the  geometry.  So,  the 
inequalities  of  (15),  with  the  help  of  the  first  inequality  of  lemma  2.4  become  (up  to  a  factor 
2  sometimes)  : 


/iP-" 


|[A.]| 


E 

caTd(I)=p 


<  h'’-"  IolI  +  eC', 


for  a  given  constant  C.  Now,  summing  up  these  inequalities  for  i  j  =  n,  with  the  help  of 
the  second  inequality  of  lemma  2.4  one  gets  our  result  for  e  small  enough.  • 

This  result  enables  us  to  detect  the  regions  of  smoothness  from  those  where  a  jump  in 
one  of  the  derivatives  occurs.  It  is  true  when  conjecture  2.1  is  true,  and  at  least  for  degree 
1  and  2. 
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3  An  E.N.O.  Reconstruction  Technique 


In  recent  papers,  Harten  and  several  other  authors  [3,  4,  5,  6,  7]  have  tried  to  derive  numerical 
methods  that  are  able  to  achieve  a  higher  order  of  accuracy  than  classical  TVD  methods. 
There  are  several  versions  of  these  techniques,  but  they  can  be  generally  viewed  in  the 
following  way  :  starting  from  some  approximation  of  a  real  function  u  (point  values  or 
average  values  in  some  control  volumes),  find  a  pointwise  high  order  approximation  u.  Two 
tools  are  then  used  : 

•  an  essentially  non  oscillatory  Lagrange  type  interpolation  of  a  function  iw, 

•  this  function  w  may  be  u  itself  if  one  starts  from  point  values  or  either  the  primitive 
function  of  u  or  its  convolution  product  with  the  characteristic  function  of  a  copy  of 
the  control  volume  if  one  can  pass  from  one  to  another  by  a  translation. 

The  latter  deconvolution  technique  can  only  be  applied,  at  least  in  its  standard  version,  to 
regulaj  meshes  as  shown  in  section  3.2. 

In  this  section,  we  want  to  adapt  both  tools  to  the  situation  of  unstructured  meshes.  At 
least  for  the  second  point,  the  situation  seems  at  first  glance  very  bad  :  the  reconstruction 
via  primitive  function  cannot  be  applied  in  the  case  of  unstructured  meshes  because  the  only 
solution  would  be  to  apply  it  to  integrals  over  domains  like  2>a/  =  [ai,6i]  x  [0.21  h]’,  possibly 
with  a  suitable  transformation  of  the  plane  as  in  [8],  which  are  not  in  general  the  union  of 
control  volumes. 

Now,  the  reconstruction  via  deconvolution  technique  can  only  be  applied  to  regular 
meshes  (i.e.  meshes  where  the  control  volumes  are  translated  form  one  node  to  another).  In 
this  section,  we  show  how  to  adapt  this  technique  for  irregular  meshes. 

In  what  follows,  T  is  a  triangulation  of  fi,  a  domain  of  IR^,  u  is  a  function  defined  on 
that  domain.  Around  each  node  i  of  T,  we  may  define  a  control  volume  C,  in  many  different 
ways.  An  example  is  (see  Figure  2)  the  control  volume  whose  boundary  is  the  segment 
joining  the  centroids  G  of  the  triangles  {i,j,k)  having  i  as  a  vertex  and  the  middles  Ii,l2 
of  the  segments  of  those  triangles  (type  I).  Another  type  of  control  volume  is  obtained  by 
considering  each  triangle  as  the  control  volume  of  its  centroid.  The  triangulation  we  need 
to  describe  the  ENO  algorithm  is  not  T  but  another  one  built  from  the  centroids  of  the 
triangles  of  T.  These  control  volumes  always  satisfy  : 

•  Ui€T  C*!  is  the  numerical  domain  Qh,  an  approximation  of  the  physical  one, 

•  Cif)Cj  is  of  empty  interior  when  i  /  j  (generally  speaking,  a  collection  of  segments  or 
an  empty  set). 
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3.1  Non  oscillatory  interpolation 

Let  us  consider  n  >  0.  In  this  section,  we  show  how  to  generalize  the  n  +  1-th  order  E.N.O. 
interpolation  technique  exposed  in  [4,  5]  to  unstructured  meshes.  As  the  results  of  chapter 
2  indicate,  we  must  deal  with  meshes  where  the  stencils  we  need  satisfy  the  property  of 
lemmas  2.1  and  2.2.  This  will  be  the  case  for  most  meshes. 

Let  T  be  a  triangulation  of  fl,  a  domain  in  IR^,  and  u  a  function  defined  on  that  domain. 
The  results  we  have  obtained  in  chapter  2  can  be  summarized  as  follows  : 

•  if  is  an  admissible  stencil  such  that  u  is  smooth  in  its  convex  hull,  then 

t+j=n 

remains  finite, 

•  if  in  the  convex  hull  of  5^"^,  u  admits  continuously  differentiable  derivatives  only  up 
to  the  order  k  <  n,  then 

Y,  i\  =  o(l“‘”>l 

«+j=n 

where  h  is  the  diameter  of  . 

Then,  as  suggested  by  [4,  5],  the  E.N.O.  algorithm  we  propose  is  to  consider  ni(u)  defined 
on  Ci  by  the  following  recursive  algorithm  : 

For  a  node  i, 

1.  Let  {Ti}  be  the  set  of  triangles  of  T  having  i  as  a  vertex.  Consider  all  the  linear 
interpolations  where  the  Tj’s  are  the  stencils.  Choose  the  one,  T^in,  where  the  sum 

»+j=i 

is  minimal.  We  set  =  the  nodes  of  T^nin- 

2.  Let  be  the  stencils  defined  at  the  previous  step.  Consider  all  the  nodes  sur¬ 
rounding  in  T  and  consider  all  the  stencils  obtained  form  by  adding  n-t- 1 

of  the  nodes  surrounding  Choose  the  stencil  minimizing  ; 

i+i=n 

We  have  intentionally  left  the  second  point  imprecise  because  it  is  obvious  that  the 
number  of  stencils  to  consider  is  in  general  huge.  To  give  an  example,  if  one  considers  Figure 
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3  for  which  n  =  2,  one  sees  that  possible  stencils  are  the  vertices  of  triangles  Tmin  and  3  of 
the  ten  other  triangles.  This  can  be  repeated  for  each  of  the  three  edges  of  Tmin  and  leads 
to  a  total  of  3  X  Cfo  =  360  possible  stencils  !  So,  one  has  to  define  criteria  for  choosing  the 
“good”  and  “bad”  stencils.  These  criteria  are  essentially  heuristic  and  a  priori  ones  . 

One  that  seems  natural  is  that  when  one  considers  the  control  volume  around  each  node, 
the  collection  of  the  control  volumes  of  all  points  of  the  stencil  should  be  convex.  Another 
one  is  that  the  criteria  leads  to  the  smallest  possible  number  of  stencils,  but  the  stencils  must 
not  be  confined  in  a  particular  angular  area  of  the  plane,  in  order  not  to  favor  any  direction. 

With  this  in  mind,  two  possible  sets  of  stencils  for  third  order  interpolation  are: 

•  the  nodes  of  triangle  Tmin  plus,  for  each  of  its  edges,  the  three  additional  nodes  of 
triangles  Ti,  Tj,  T^.  This  leads  to  a  maximum  of  three  stencils  per  triangle, 

•  or  the  nodes  of  triangles  of  Tmin  plus,  for  each  of  its  edges,  the  three  additional  nodes 
of  triangles  of 

-  ri,T2,T3, 

-  Ti^Ti  and  T4  or  T5, 

-  Ti,  Ta  and  Te  or  TV, 

-  Ti,  Tg  or  Tio  and  one  of  the  six  triangles  Tz,  T3,r4,  Ts^Te^Tr. 

The  second  solution  leads  to  a  maximum  number  of  52  stencils  once  Tmin  has  been  found. 
We  have  made  several  tests  to  evaluate  the  “performance”  of  each  type  of  stencil.  They  are 
given  in  section  4. 

This  particular  interpolation  is  n  +  1-th  order  accurate;  because  u  is  a  polynomial  of 
degree  less  than  n,  we  have  ni(w)  =  u.  This  property  ensures  the  n  -f  1-th  order  accuracy 
[15],  and  in  particular,  we  have  the  estimations  of  theorem  2.3. 


3.2  Deconvolution  technique  revisited 

If  (ar),et^  is  a  regular  mesh  of  R  and  u  is  a  real  valued  function  on  IR,  the  reconstruction  by 
the  deconvolution  technique  consists  of  applying  the  previous  algorithm,  not  to  u  but  to 

Ux{y)  = -  /  u{x-Vy -Xi)dx,  (16) 

where,  as  usual,  the  mesh  size  Ai  =  Xi+i  —  Xi  is  constant  and  ii4.i/2  =  x,  -f  Ax/2.  In 
particular,  we  see  that  Ui  does  not  depend  on  i  and  that  u{xi)  is  the  average  of  u  on 
a:,+i/2]-  These  values  are  assumed  to  be  known.  Let  ni(u)  be  the  m  -f  1-th  order 
Lagrange  interpolation  as  decribed  in  the  previous  section,  with  m  >  n.  Then,  the  idea  is  to 


perform  a  Taylor  expansion  of  u  and  its  successive  derivatives  around  x,-,  to  truncate  them 
at  order  n  —  A:,  to  replace  the  values  u,  •  •  •,  by  ni(u)  and  its  n  successive  derivatives, 
and  to  replace  the  values  of  u  by  those  of  n2(tt),  the  approximation  of  u  we  are  looking  for  : 


where 


I  dx 

-  •’^,-1/2 _ 


The  linear  system  is  easily  invertible  because  the  matrix  is  upper  triangular  and  its  diagonal 
consists  of  all  “l”s.  Furthermore,  it  is  shown  in  [4],  for  example,  that  the  average  value  of 
n2(u)  over  [x,_i/2,  x,+i/2]  is  exactly  Uj.  Last,  this  approximation  has  the  desired  order  of 
accuracy  when  u  is  smooth  because  polynomials  are  left  invariant  by  the  construction. 

This  latter  point  is  the  fundamental  reason  why  one  achieves  the  expected  order  of 
accuracy  .  Polynomials  are  left  invariant  by  the  construction  because  the  shape  of  the 
control  cells  does  not  change  from  one  point  to  another.  If  this  where  not  the  case,  i.e  if 
^^i+i/2  —  ^i+i  ~  were  not  constant,  the  formula  (16)  would  indeed  depend  on  the  point 
Xi  and  this  property  would  be  lost.  In  order  to  show  this,  we  simply  consider  u{x)  =  x  and 
a  m  +  1-th  order  interpolation  that  has  values  u,  at  points  x,.  Assuming  that  {xq,  Xi,  •  •  •}  is 
the  stencil  selected  by  the  E.N.O.  algorithm,  we  have  : 


ni(u)  =  «o  +  Kix  -  Xq)  + 


where  ^  ^ 

2  Axi/2 


When  the  mesh  is  not  regular,  K  ^  I  m  general.  To  obtain  n2(n)  =  n,  one  must  have  : 

ni(u)(y)  =  u{y)  -b  aiu'(y), 
n2(w)(y)  =  K  =  u'{y). 

The  second  equation  indicates  that  one  must  have  K  =  I  which  is,  in  general,  not  true. 

To  overcome  this  problem,  we  propose  the  following  technique  :  apply  the  ENO  search 
algorithm  not  to  u  defined  by  equation  16  but  to  : 


u{y)  = 


area{Cg( 


•'C  r„) 


X  -f  y  -  Xo)dx, 
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where  5^"^  is  any  possible  stencil  around  the  node  xq  that  one  has  to  test  and  Cs(n)  is  the 
union  of  the  control  volumes  of  each  node  in  : 

U  Ci-  (19) 

jesf.") 

Now,  one  has  to  evaluate  the  integral  (18)  from  the  average  value  of  u.  This  cannot  generally 
be  achieved  for  any  function,  but  is  possible  for  the  polynomials  of  IRn[X,  F],  at  least  in 
general.  This  will  be  true  for  the  N(n)  linear  forms  over  IRnfA',  F]  because 

<P>i=-^--l  P{x)dx,  forall  PelRn[X,F],  (20) 

are  independent.  For  all  the  meshes  we  have  considered,  these  linear  forms  where  always 
independent,  so  the  problem  had  a  solution.  If  this  is  true,  then  one  can  find  coefficients 
1  ^  f  ^  N(n)  so  that 

1  f 

- — - -  /  u{x  +  y-  xo)  dy=Yl  ^  (21) 

area[Gg{n))  /_j 

when  u  belongs  to  IR„[X,  F].  If  not,  the  equation  (21)  gives  a  n  +  1-th  order  quadrature 
formula.  With  all  this,  we  get  the  following  theorem  whose  proof  is  obvious. 

Theorem  3.1  The  algorithm  defined  by  the  E.N.O.  technique  with  (18),  (19),  (20)  and  (21) 
leaves  invariant  the  polynomials  o/IRn[X,  F]  and  hence  gives  an  +  1-th  order  approximation 
of  smooth  functions.  Moreover,  this  approximation  is  conservative  up  to  the  quadrature 
formula  21. 

3.3  Some  remarks  for  the  practical  calculation  of  the  reconstruc¬ 
tion 

In  section  2.1.3,  we  have  discussed  the  problem  of  the  practical  determination  of  the  coeffi¬ 
cients  of  a  Lagrange  interpolant  because  the  linear  systems  to  be  solved  have  in  general  small 
coefficients.  The  same  problem  also  arises  here  if  one  uses  the  monomials  (X— Xo)’(F — yo)^  to 
determine  the  a/(y)  of  equation  21.  This  can  easily  be  seen  by  using  Hadamard’s  inequality 
as  in  section  2.1.3. 

For  a  given  node  xq,  the  E.N.O.  technique  we  propose  naturally  introduces  one  triangle 
having  that  node  eis  vertex,  the  triangle  Tmin  as  in  Figure  3  .  So,  as  in  section  2.1.3,  we  will 
use  the  barycentric  coordinates  towards  that  triangle  for  practical  computations. 
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4  Numerical  examples 


We  have  performed  several  tests  on  the  second  and  third  order  E.N.O.  interpolation  and 
E.N.O.  reconstruction,  but  we  only  report  the  third  order  results  since  they  are  a  priori 
more  challenging.  In  particular,  we  intend  to  check  numerically  that  the  expected  order  of 
accuracy  is  in  fact  reached  for  smooth  functions. 

The  two  types  of  stencils  have  been  tested  on  the  latter  case.  The  use  of  the  second  type 
of  possible  stencils  results  in  a  much  more  expensive  approximation  (in  general,  one  must 
test  52  stencils  per  triangle  versus  only  3  in  the  simplest  version)  and  the  results  have  never 
been  dramatically  improved.  All  the  results  that  are  presented  bellow  have  been  obtained 
with  the  3  stencil  version  of  the  method. 

In  all  these  examples,  we  have  assumed  that  the  control  volumes  are  of  “type  I”.  The 
practical  calculations  of  the  averages  in  these  control  volumes  have  been  performed  with  a 
5-th  order  quadrature  formula  [14]. 

The  tests  on  smooth  functions  will  be  performed  on  : 

u{x,y)  =  cos{2t{x^  -\-  y^)). 

We  have  displayed  in  Figure  4  the  error  of  the  interpolation.  The  plain  curve  with 
squares  is  obtained  with  the  E.N.O.  interpolation,  the  plain  curve  with  circles  is  obtained 
with  the  E.N.O.  reconstruction.  The  dashed  line  indicates  the  slope  —3.  One  can  see  that 
the  expected  order  of  accuracy  is  indeed  reached.  The  mesh  size  h  has  been  measured  by 
choosing  the  largest  segment  of  the  triangulation.  All  the  error  estimates  have  been  obtained 
on  irregular  meshes  as  the  one  presented  on  Figure  5.  These  meshes  are  obtained  as  random 
perturbations  of  regular  structured  meshes.  The  set  of  points  that  one  obtains  is  triangulated 
by  the  Brower  algorithm  to  get  a  Delaunay  triangulation.  The  main  difference  between  such 
a  mesh  and  the  regular  structured  one  is  that  the  number  of  triangles  each  node  belongs  to 
is  different.  We  also  have  done  the  same  tests  with  regular  meshes,  and  we  have  not  seen 
any  degradation  of  the  convergence. 

The  locally  smooth  function  we  have  chosen  is  obtained  by  a  modification  of  that  used 
by  Harten  in  [7]  for  example  :  if  is  any  angle,  let  be  : 

I  if  r  <  U{x,y)  =  -rsin(|r2), 

y)  =  ■!  if  r  >  |,  /^(x,  y)  =  2r  -  I  sin  (37rr),  where  r  =  i  -|-  tan  {<f>)y,  (22) 

I  if  H  <  h  =  |sin(27rr)l, 

and  let  u  be  : 

{if  x<^cosiry,  u{x,y)  =  f^^{x,y), 

if  x>2COS7ry,  u(i,  y)  = /_^y^(x,  y) -f- cos  (27ry). 
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(23) 


The  function  defined  by  (22)-(23)  shows  discontinuities  in  the  function  itself  and  its  first  order 
derivatives;  some  of  the  discontinuities  are  straight  lines  (never  aligned  to  the  mesh),  one  is 
a  curved  line  where  the  jump  changes  from  one  point  to  another.  Last,  the  behaviour  of  u  is 
basically  one-dimensional  on  the  left  of  the  curve  x  =  cos7r7//2  and  really  two-dimensional 
on  the  right. 

A  plot  of  this  function  is  given  in  Figure  6.  One  should  obtain  straight  lines  and  smooth 
discontinuity  transitions  contrary  to  what  is  shown  in  the  Figure  :  this  is  an  effect  of  the 
graphic  device  adapted  to  Pi  interpolation. 

In  Figures  7  and  8,  we  have  displayed  the  node  values  of  the  E.N.O.  reconstruction  for 
two  meshes  (1600  nodes  and  6400  nodes).  To  better  see  the  behaviour  of  both  approximation 
techniques,  we  also  present  cross-section  on  three  lines  :  Y  =  0.75,  F  =  0  and  Y  =  —0.45. 
The  approximations  are  obtained  from  the  1600  nodes  mesh  (see  Figure  5).  The  latter 
line  goes  through  one  of  the  triple  points  (see  Figure  6).  One  can  see  that  the  various 
discontinuities  and  the  smooth  regions  are  well  captured  by  both  techniques. 

To  end  this  section,  we  must  note  that  the  algorithm  for  choosing  the  stencil  may  lead 
to  some  difficulties  at  the  boundaries  as  can  be  seen  in  Figure  6  on  the  left  upper  corner  : 
the  most  left  upper  triangle  of  the  mesh  (Figure  5)  does  not  admit  any  additional  points  of 
the  type  we  consider  to  make  a  stencil. 

5  Conclusion 

In  this  report,  we  have  developed  two  methods  for  the  reconstruction  of  a  function  admitting 
discontinuities  only  on  regular  planar  curves  from  their  node  values  or  from  their  average  in 
control  volumes  that  surround  them.  In  order  to  give  a  firm  basis  to  the  Essentially  Non- 
Oscillatory  interpolation,  we  have  studied  the  behaviour  of  the  highest  order  coefficients  of 
the  Lagrange  interpolation  of  smooth  functions  and  unsmooth  ones  for  which  the  discontinu¬ 
ities  lie  on  regular  curves.  We  have  also  given  an  adaptation  of  the  so  called  “reconstruction 
via  deconvolution”  to  irregular  triangulated  meshes. 

These  techniques  have  been  shown  to  work  quite  well  on  smooth  and  unsmooth  functions. 
In  particular,  we  have  shown  in  these  examples  that  the  minimum  number  of  possible  stencils 
was  sufficient  for  our  purpose. 
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Figure  2:  Control  volume  around  i 
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Figure  3:  Some  possible  interpolation  points.  Circles  :  points  of  Tmin  (second  order  in¬ 
terpolation),  black  circles  :  points  that  may  be  added  to  obtain  a  stencil  for  third  order 
interpolation. 
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Figure  4;  L°°  error  for  /(x,  y)  =  cos  [27r(x^  +  y*)].  Squares  ;  E.N.O.  interpolation  only, 
Circles  ;  E.N.O.  +  reconstruction.  Dashed  line  :  slope  -3 
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Figure  5:  Typical  mesh.  1600  nodes,  3042  triangles. 
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Figure  6:  Exact  function,  Mesh  with  6400  nodes.  Min=-1.331,  Max=2.650,  6  =  8.292  10"^. 
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Figure  7:  E.N.O.  reconstruction  with  the  1600  nodes  mesh.  Min=-1.308,  Max=2.651,  6  = 
8.249  10-2. 
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Figure  8;  E.N.O.  reconstruction  with  the  6400  nodes  mesh.  Min—1.325,  Max-2.650,  8 
0.8281  10"^ 
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Figure  10:  Cross-section  at  K  =  0.75  of  the  E.N.O  interpolation  (a)  and  the  E.N.O.  recon 
struction  (b)  for  the  1600  nodes  mesh. 
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